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A Simulation to Study Speed Distributions in a Solar Plasma 


ABSTRACT 

We < Peter Cheeseman of NASA Ames/Caelum Research) & Jose Luis Alvare Uos of the 
SJSU Physics Department/SJSU Foundation.) have carried out a numerical simulation of a 
plasma with characteristics similar to those found in the core of the Sun. Particular emphas.s 
is placed on the Coulomb interaction between the ions and electrons, which could result in a 
relative velocity distribution different from the Maxwell-Boltzmann (MB) distribution 
generally assumed for a plasma. The fact that the distribution may not exactly follow the MB 
distribution could have very important consequences fora variety of problems in solar 
physics, especially the neutrino problem. Very briefly, the neutrino problem is that the 
observ ed neutrino detections from the Sun are smaller than w hat the standard solar theory 
predicts [1 ]. In Section 1 we introduce the problem and in section II we discuss the approach 
to try to solve the problem: i.e., a molecular dynamics approach. In section 111 we provide 
details about the integration method, and any simplifications that can be applied to the 
problem. In section IV (the core of this report) we state our results, first for the specific case 
of 1000 particles and then for other cases with different number of particles. In section V we 
summarize our findings and state our conclusions. Sections VI VII and Vlll provide the list 
of figures, reference material and acknow ledgments respectively. 

I. introduction/statement of the problem 

A plasma is a hot. ionized gas which is composed of (positive) ions and free electrons 
distributed over a region of space (Chen (2]: Spitzer [3]). Plasmas are electrically neutral 
(overall charged). Examples of plasmas are the ionosphere, which is an upper layer of the 
atmosphere, and the gas in the interior of the Sun. In the core of the Sun, temperatures are so 
high that atoms are completely stripped of their electrons, thereby constituting a plasma, 
composed of (positive) nuclei and free electrons. 

We wish to study the behavior of a plasma under the same conditions as in the Sun s 
core. The main interaction affecting the behavior of the ions and electrons in plasma is the 
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Coulomb Force. Much of the work done in Plastna Physics/Solar Physics assumes a MB speed 
distribution for the ions and electrons. The MB speed distribution is 

g(v) = 4rt,T(m/2;z£ fl r)' : v : exp(-v : /v" ) (1 ) 

where n is the number density, kg is Boltzmann's constant. T the temperature and we 
define v m = 2k H T ! m . The quantity g(v)dv gives the number of particles per unit volume with 

a speed between v and v + dv. Analytical derivations of the MB distribution assume that all of 
the energy is kinetic energy, that is, the classical derivation neglects any inter-particle forces 
such as the Coulomb force. A measure of the deviation from “ideal gas" conditions is the 
parameter r which will be defined later: suffice it to say that for an ideal gas. /"= 0. and for a 
real plasma f> 0. There will be deviations from MB due to the Coulomb interaction, the 
question is how much (Swihart [4]). For another view on the effects of a non-MB distribution 
on solar neutrino rates see Clayton (refs. [5] and [6]). 

Energy production in the Sun is by thermonuclear reactions, primarily via the proton- 
proton (PP) chain. There are three possible PP chains (PPI, PPII and PPIII) and which 
dominates depends primarily on temperature (Clayton. 1983, [7]). In the overall view, the PP 
chain amounts to the combination of four protons to produce a 4 He nucleus, tw o positrons and 
two (electron) neutrinos (Bahcall 1989 [8]): 

4'H — > 4 He + 2e~ + 2v e [Q = -26.7 Me V] (2) 

Of the 26.7 MeV released, only about 0.6 MeV is carried away by the neutrinos. The reaction 
is not as straightforward as implied by the above equation, which is a simplified view of the 
whole chain. The PP chain involves several "branches", such as the first branch of the PP 
chain: 

'H + 'H-> : H + e + + v e [Q = 0.42 MeV] (3) 

For the rest of the branches see (Clayton, 1983 [7]) or (Bahcall, 1989, [8]). 





The alternative reaction we are interested in is (Bahcall & Wolf. 1963 [9]) 


3 He + e' — ► 3 H + v e [Q = -18.6keV] (4) 

This is an electron capture by a 3 He nucleus to give a Tritium nucleus plus a neutrino. This 
reaction is not part of the standard solar model: note that it is endothermic and w ould release 
neutrinos of very low energy, lower than present "neutrino telescopes" would detect. If this 
reaction were to "go", the energy would have to be provided by the electrons, w hich would 
need an average energy of 18.6 keV. At the core of the Sun (1.56xl0 7 K) the mean kinetic 
energy of electrons is 2 keV. Assuming a MB distribution the ratio of electrons w ith energy 
18.6 keV to those w ith 2 keV is approximately 1.3 lxl O' 5 . In Ref. [9], Bahcall & Wolf imply 
that the reaction ? He(e'.v e ) ; ’H is only important in the case that the central density of a star p it > 
2x 1 0 7 Kg/m 3 , and so would proceed at a negligible rate in the Sun's core, where p = 1 .48x 1 O' 
Kg/m 3 . However, they did not take into account the statistical mechanics of interacting 
particles. We suspected that the mean energy of electrons may be higher than 2 keV, w hich 
would make the electron capture reaction given above to be much more likely than the 
standard solar model predicts. A possible cause for deviations from MB distribution are the 
Coulomb interactions between electrons and ions. Another possible cause is electron 
degeneracy (see Swihart, 1972 [4]), but we concentrate on the effect of Coulomb interactions 
here. This is the purpose of the simulation. 

II. APPROACH TO THE PROBLEM 

We have dev eloped a numerical model of a plasma, with characteristics ty pical for the 
interior of the Sun; i.e., the plasma w ill have the same temperature and density ( T - 1 .56x10 
K. p- 1 .48x1 0 5 Kg/m 3 , see Bahcall [8]), as the core of the Sun. We used a Molecular 
Dynamics (MD) simulation of electrons and ions (Haile [10]), using paired electron-ion 
"clusters" (as dipoles) to cut off field contributions below a given threshold. We use the 
Velocity-Verlet scheme to advance the system in time. It is reasonably fast and quite accurate. 
The simulation consists of a large number of ions and electrons in a given volume (determined 
by the density in the Sun’s core) using periodic boundary conditions. We ignore quantum 
effects. We are investigating the following (among other things): 



1 . The possibility of a non-Maxwellian velocity distributions (ions and electrons) due to 
the long-range effects of the Coulomb force. This (possibly) non-Maxwellian 
distribution can affect the rate of certain nuclear reactions occurring in the Sun s core 
which are very important for the solar neutrino problem. 

2. Check the relaxation times for ion-ion and ion-electrons interactions. 

3. We will also check the (analytical) mean-free-path formula for ions and electrons. 

III. DETAILS OF THE INTEGRATION 
(a) Units 

We choose to use atomic units or au (see McQuarrie [1 1]) as the most natural set of units 
for the problem at hand. The unit of length is the Bohr radius, the unit of mass is the electron s 
mass. etc. See the following table: 


Table I: Atomic Units 


Quantity 

Atomic Unit 

MKS Equivalent 

Mass 

m= 1 (electron mass) 

9.1091x10 Jl Kg 

Charge 

|c|=l (electron charge) 

1.6021xl0' lv C 

Length 

a a = 1 (bohr) 

5.2918x10 " m 

Time 

Electron period in Bohr orbit 

2.4189x10 *' s 

Energy 

1 hartree 

4.35944x10 '* J 


Only temperature is the same as in the MKS set, i.e.. it is measured in degrees Kelvin. 

(b) Initial Conditions 

Since one of the things we w ant to check are the relaxation times, we choose initial 
conditions which are as far removed from equilibrium as possible, and then we see how long it 
takes for the system to relax. The initial positions are arranged in an FCC (Face-Centered- 
Cubic) pattern, as in a NaCl crystal, with electrons and protons at alternate positions. The total 
number of particles N p is constrained by the fact that we want the overall charge to be zero. If 
we assume we have only electrons and protons, as in our case, this means that the total number 
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of particles has to be an even number. Furthermore, since we are using an FCC pattern, the 
cubic root of the total number of particles has to be an integer: this is the number of particles 
on a side of the computational box. Therefore. N p is constrained by the following conditions: 

• X p is even 

• (X p )' J is an integer 

Using these criteria we find that N p — 8/ 3 , with /— 1 .2.3. 4. 5. etc., so Ap = 8, 64. 216.512. 1 000. 
1728. etc. 

The initial spacing d between each electron and proton in the lattice is determined from 
the Wigner-Seitz radius a , or 



where n, is the ion number density (equal to the electron number density ). If we assume only 
Hydrogen is present, w e obtain a = 0.26. As prev iously mentioned, a measure of the Coulomb 
coupling strength is T, which is defined as the ratio of the average potential energy to the 
average kinetic energy, or for equilibrium conditions 


f = 


2 

3kja 


( 6 ) 


Using the previously obtained value of a. w e find 7~= 0.05 (in au ke~ 3. 17606x1 0" 6 
hartrees/K). 

Since w e are using a cubic computational box. it can be shown that the length of a side 
will be L = dN ^ , so that the volume is V = d'N p . The number density of ions n, is 


number of protons \ 

Volume of comp, box d'N f 2d 5 


( 7 ) 


Substituting this expression into Eq. (5) and solving for d we obtain 
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a 


( 8 ) 



Once we ha\e decided the total number of particles and their initial positions we need to decide 
on their initial velocity. At first thought this seems simple: assign the initial velocities 
according to the average kinetic energy per particle at the center of the Sun (determined by the 
temperature). However, the average potential energy per ion-electron pair is lower at 1=0.0 
(for an FCC lattice) than for equilibrium conditions. This means that as the system evolves in 
time and approaches equilibrium, its potential energy will increase slightly. To conserve 
energy, the kinetic component has to decrease as the system evolv es, which means the initial 
kinetic energy has to be slightly greater than the kinetic energy at equilibrium. We want the 
total energy per particle (after relaxation) to be equal to the mean kinetic energy per particle at 
the center of the Sun. The extra amount of kinetic energy is obtained by comparing the 
potential energy per ( proton-electron ) pair at t- 0.0, u(0), and at equilibrium. <u>. At / = 0.0 
we have an FCC lattice, and the potential energy per pair is (in au’s) 


( a ) 

m 

v« J 

V2 n) 


X 


S -1.3659 fa 


(9) 


where a= 1.7476 is the Madelung constant for an FCC lattice (see Ashcroft & Mermin [12]). 
At equilibrium, the average distance between an ion and an electron is the Wigner- Seitz radius, 
so 

<u>= -Ma (10) 


The initial kinetic energy per particle k(0) w ill be the kinetic energy' at equilibrium <k> plus 
the difference between Eq. (9) & (10) divided by 2: 


k( 0) =< k > + 


< u > -w(0) 
2 


( 11 ) 
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where < k >= 3k „T ! 2 . Since all particles have the same initial kinetic energy, the initial 
speed per particle will be v(0) = -J2k{0)/m . Electrons will be faster than protons (by a factor 

of - 43) due to their lower mass. In effect, the initial speed distribution for electrons and 
protons is a delta function. The direction of the velocity vectors is chosen at random. 

(c) Advancing the System in Time 

Once the initial conditions are chosen (in the form of initial positions and velocities), the 
system is advanced in time via the Velocity-Verlet scheme (Allen & Tildesley 1992 [13] or 
Swope et al. [14]). Each particle /= 1 ,2,3 \ p is advanced in time as follows: 

r.tl-At) - r,(i) + v,(t) »Al - 0.5*[F,(t) m,J*Ar ( 1 2a) 

v , (t-At 2) = v,(t) - 0.5 •[F/t) m,J»At ( 1 2b) 

v, (t~At) = v, (t+At/2) - 0. 5 •[F,(t+Al). mj *J/ (12c) 

where At is the (fixed) time-step size and (F/mJ is the acceleration on particle /'. The scheme 
is reasonably fast and very accurate: the local truncation error is 0(4/"'). We use periodic 
boundary conditions: once a particle exits the computational box. it re-enters it on the opposite 
side. 

(d) Forces 

The computation of the forces on particle i is done as follow s. A sphere centered on 
particle / with a ‘‘cut-off radius R c is chosen so that the forces on particle / of (most) particles 
outside of this sphere are ignored. Choosing the quantity R c is a compromise between 
execution speed and accuracy, and must be found empirically; we chose to use Rc = 0.67. In 
addition, we keep charge neutrality for all the particles whose forces are acting on particle ; 
(including particle i ). If there is no charge neutrality within this sphere, we search outside the 
sphere for the nearest particles of the appropriate charge to keep the net charge zero. 
Thereafter, the force on particle / is 
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(13) 


j m * 


where the sum goes from j =1.2.3 V„. where .V„ is the number of "neighbors" of particle / 

(which varies as a function of time, but typically <,V„> - 40 for Rc - 0.67) and F-,j is the force 
on particle / due to particle j. The reason for the emphasis on charge neutrality is that we want 
to use paired electron-ion "clusters" (as dipoles) to cut off field contributions below a given 
threshold, as the field of a dipole falls off faster than that of a simple point charge. Finding 
w hich particles to use in the calculation of the forces on a particle (the "neighbors") is 
computationally expensive, so we use the concept of neighbor lists to save time (see Haile, 
[10]). A list of the N„ neighbors of particle / is obtained and the same list is used over several 
consecutive time-steps to compute the forces. Typically we update the list even 15 time-steps 

We use a modified form of the Coulomb interaction between two particles: 


f = 


MS* 


( r ';+£ z f 


(14) 


which, in the limit of small e tends to the real Coulomb force. The main reason for this 
"softened" force is one of computational convenience, as it removes the possibility of the 
denominator in Eq. (14) from "blow ing up" due to very close approaches. Another way to 
think of this softening is that each particle has a cloud of charge, 85% of which is contained 
within a radius of3f (see Sellwood [15]). Again, one must choose empirically a small enough 
value of £so that Eq. (14) is a good approximation to the real Coulomb force, but not so small 
that close approaches between particles result in large errors (because of the finite size of the 
time-step At) due to extreme values of the denominator. We use e- 2x10°. 

The force between a "real" particle / (located inside the computational box, CB) and an 
“image” particle j’ (a projection of a “real” particle j close to the opposite edge of the 
computational box from particle i) is handled by Haile’s “ Minimum Image Criterion' (see 
[10]). This becomes important for particles close to the edge of the CB. This criterion 
identifies which of the multiple images of a real particle will be used in the computation of the 
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forces, and is applied separate!) to each cartesian component. For example, for the x 
component we have 


x^x ; ,-L (x, <-L/2) 

( x , > L / 2 ) 


and x y is unchanged if - 1/2 < x r < 1/2 . 

Finally, we mention that another computational “trick" to speed up the simulation is to 
use Newton's Third Law whenever possible. If the force on particle j due to particle / F p is 

computed, then the force on particle i due to j F u does not need to be computed again as it is 
already known, i.e., F n = -F jt (along the same line!). 


(d) Stochastic Heating 

Stochastic Heating is a numerical artifact and is a result of all the errors in the 
computer simulation (due to truncation error, finite time-step, rounding off, finite computer 
precision, etc.). The result is that the mean kinetic energy per particle <k> increases as the 
simulation proceeds in time (see Hockney & Eastwood [16]). The average error in the kinetic 
energy' <h> of, say. particle /. is a linear function of the number of time-steps n. and is given 
by 


<h>= 


2m, 


05) 


where q, is the particle's charge, m, is its mass, S is the (modulus of the) error in the electric 
field (due to truncation, finite machine precision, etc.) computed at the location of particle / 
and At is the time-step used. It can be seen from Eq. (15) that the smaller the time-step, the 
smaller the effects of stochastic heating. Furthermore, the larger the mass, the less error in 
kinetic energy; hence, for a given time-step size, protons are much less affected (by a factor of 
m p™anl m 'kcmm ~ 1836) than electrons by stochastic heating. 
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The result is that unless we pawide corrections, the kinetic energy per particle (and 
hence the total system energy) will increase without bound. Therefore, we periodically check 
the error in the total energy and automatically provide corrections for stochastic heating by 
reducing the magnitude of the velocity vectors of all the particles. These corrections are 
different for protons than electrons. 

(e) Table of characteristic lengths 

Follow ing is a table of the characteristic length scales involved in the problem. Also 
included is the Debye length, often-used in Plasma Physics work to indicate the "screening" 
effects of a collection of charges: 


Table II: Characteristic Lengths for the problem 


Length 

Value (au) 

Wigner-Seitz radius, a. Eq. (5) 

0.26 

Lattice dist., d. Eq. (8) 

0.33 

— | 

Cut-off radius. R c 

0.60 <R C < 0.67 

Softening radius. £ 

2.0x1 O' 3 

Debye length. Ad 

0.54 


The significance of the Debye length is that for distances much greater than /.£> from a giv en 
charge, screening effectively cancels out the charge (see Chen. [2]). In addition, if the number 
of particles inside a "Debye" cube is much greater than unity, then the sy stem is defined as a 
"collision-less" sy stem (see [16]). In our case, the total number of particles in a Deby e cube is 

( X D jdt) > * 4 . Therefore, collisions are important for our problem, as expected. 

(f) Running the Programs 

We implemented all of the previously mentioned ideas in the form of FORTRAN 77 
codes. There are actually two v ersions of the program. The first version (main: 
PLSMMD.FOR) advances the system in time from t = 0.0 (FCC lattice, initial speed 
distribution is a delta function) to a sufficiently long enough time that the sy stem is "well 
relaxed”. Along the way the program checks the speed distribution, conservation of energy 
and conservation of momentum. A time after which the system is “well relaxed" has been 




found empirically and it is approximately 330.0 time units (= 330.0 x 2.4 189x1 O’ 1 7 - 8x1 O' 15 
seconds). The second version (main: PLSMMD2.FOR) integrates the system in time from 
where the first program left off to an additional, arbitrary number of time steps. In this second 
version we check the speed distribution, energy conservation and close approaches. To speed 
up the simulation we skip momentum checking. The simulations are quite CPU intensive and 
a typical run could have lasted anyw here from a few days to a month depending on the total 
number of particles. The computer used w as a 6-CPU. Ultra Enterprise 4000/5000 Sun 
Micros> stems workstation w ith 1.2 GB of main memory running Solaris5.6. Both versions 
were compiled with Sun Microsystems' f77 compiler in double precision and optimized for 
speed (i.e.. used the [-04 -fast] compiler options). Each run took place on a single CPU. a 
248 MHz UltraSPARC-11 chip. See the Appendix for the programs and a discussion of the 
input parameters for these two programs. 

IV. RESULTS 

We have done several runs with different numbers of particles. In what follows we 
will discuss a typical case (with N p =1000 particles), state our results and then tabulate the 
same quantities using other values of A' p . Whenever possible we compare our answers w ith 
analytical expressions. 

(a) Case N p =1000 particles (Relaxation Run): 

We ran the first program. PLSMMD at the Unix prompt. This first program is run 
mainly to "relax" the system to equilibrium conditions, so we call this a "relaxation run". In 
Figure 1 we see a plot of the initial positions and velocities in the computational box (CB). 

The length of a side the CB is L = dN p ' \ so for \' p = 1000 we have L = 3.328. The 

arrangement of electron-proton positions can easily be seen to be that of an FCC lattice: the 
circles w ith dots at the center are protons while the dots are electrons. The velocity vectors are 
also shown, although at this scale only the electron velocity vectors can be seen (the velocity 
vectors of the protons are approximately 43 times shorter). Note that the velocity vectors are 
randomly oriented, but the magnitude of all these vectors is the same (for a given species, 
obtained from Eq. (11)). For electrons this initial speed is v w (0) = 12.2322 while for protons it 

is v(0) =0.2854. In effect, the initial speed distribution is a delta function. It's easy to see 
that, apart from the orientation of the velocity vectors, there’s nothing random about these 
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initial conditions. We now propagate this system forward in time using the Velocity -Verlet 
scheme. The time step size is At = 1 . 1 x 1 O' 4 , integrated for 3x 1 0* time-steps. The softening 
radius is e- 2x1 O' 3 and the outer "cut-ofT radius R c = 0.67. 

In Fig. (2) we see four panels with information related to the behavior of the energy for 
this first integration. In the first panel (Fig. 2(a)) w e see the behav ior as a function of time 
( 0.0 < t < 330.0 ) of the total kinetic energy K(t) (top line), the total sy stem energy, both 
analytical and numerical. £ and E„(t). (very close together appearing as one line, right below 
the kinetic energy) and the total potential energy U(t) (bottom line, always negative). Note that 
most of the energy is kinetic, and it is this fact that makes plasma researchers ignore the 
potential energy of the sy stem in some cases, assume ideal gas conditions, and hence use pure 
Maxwell-Boltzmann (MB) statistics. The two circles on the vertical axis represent the 
analytical values of the initial potential and kinetic energies (U 0 = u(0) *\' p 2 = -2626.8 
hartrees. K 0 = kf0)*N p - 74813 hartrees. with u(0) and k(0) given by Eqs. (9) and (11) 
respectively) from which we can obtain an analytical value for the total system energy £ - 
721 86.2. In the second panel, Fig. 2(b), we have zoomed in on the first panel to observe the 
behavior of the system energy at the beginning of the simulation (note that the time axis goes 
from 0.0 to 5.0). We also moved U(t) upward by 72000 hartrees to better appreciate the 
behavior of the kinetic and potential energies as the sy stem moves towards equilibrium. Firstly 
notice that the (total) kinetic and potential energies correctly keep “sync" w ith each other to 
approximately conserve the (numeric) total energy E„(t) = K(t)~ L (t) - const.. The straight 
line is £ and the (not so) straight line just below it is E„(t). Second, note how the potential 
energy moves upwards as time increases. Again this is due to the fact that at t = 0.0 we have 
an FCC lattice, which has a lower potential energy than a random configuration. 

Consequently, the potential energy tends to increase asymptotically towards an equilibrium 
value and the kinetic energy tends to decrease asymptotically to it own equilibrium value 
(corresponding to the system temperature). In the third panel Fig. 2(c), we see the proton and 
electron "components” of the kinetic and potential energies. Note that the potential energies of 
electrons and protons are very close to each other (indeed, they are almost indistinguishable at 
this scale), but the kinetic components are quite separated. The electron kinetic energies are 
higher than the proton kinetic energies, or in other words, the electrons are at a higher 
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temperature than the protons; however, note also that we started all particles with the same 
kinetic energy. In Fig. 2(d) we plot the percent error in energy, which we define as 

%error = 1 00 x EAll — — ( 1 6) 

£ 

It can be seen that the |%error; < 1 % during the simulation. In fact, the worst error was ~ 0.9 
°o at / = 93.2 and the mean of the absolute value was <|%error> - 0.22%. Also worth 
mentioning is that we have checked the numerical value of 7" and we obtained <F> = 0.0521 . 
in good agreement with the discussion given previously (see the discussion regarding Eq. (6)). 

In Fig. (3) we see the behavior of the system momentum P(t) as a function of time. 
Since there are no external forces acting on the system, the system momentum should be 
conserv ed. In Fig. 3(a) we plot the total system momentum as a function of time. We plot 
here the x, y and z components of the total momentum [p x . p y , p : ) as well as the magnitude of 
the total system momentum p lol where 

p m = ip) + p; + p: )' : ( |7 ) 

Note that conservation of momentum holds quite well: not only is the total momentum well 
conserv ed but so are each of the three components. In Fig. 3(b) we see a plot of the "electron" 
component of the momentum and Finally on Fig. 3(c) we see a plot of the "proton" component 
of the momentum. It can be seen that most of the momentum is carried by the protons (due to 
their higher mass). 

In Fig. 4 we see four panels showing data for a specific particle, in this case an 
electron. In the first panel (Fig. 4(a)) w e see a plot of the electron’s position as a function of 
time in the form of x(t), y(t), z(t). At t - 0.45 z(t) undergoes a large deflection, then y(t) 
undergoes a large deflection at / ~ 0.75 and finally z(t) also undergoes a large deflection at t- 
1 .2. Because we are using periodic boundary conditions the following constraint holds: 

)x(/)j < £/ 2 = 1 .664 where L is the length of a side of the CB; the same relationship holds for 
yftf, z(t). In Fig.4 (b) we see a plot of the velocity as a function of time in the form of 
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V r (/).V (t).v.(t) . It can be seen that at 1=0.0 , the z-component was the largest component. In 

the follow ing panel Fig. 4(c) we plot the cosine of the angle ^between this electron's initial 
velocity vector v 0 = v(t = 0.0) and its velocity vector as it evolves as a function of time v = v(t), 
i.e., the relationship between these two vectors and 0 is just v a » v = cos(#). We are interested 
in measuring the deflection time, defined as the time "in which deflections deflect the test 
particles by 90°" (Spitzer [3], p. 1 3 1 ): by definition this is the time when cos{0) = 0.0. In this 
case we found to — 0.80. Finally, we define the Mean Free Path A as the distance the particle 
travels during the deflection time: in Fig. 4(d) w e plot the particle's speed 
v(/) = [v;(/) + v*(f) + v; (/)] i: as a function of time. The integral of vft) from 0.0 to to gives us 

the electron's mean free path; here A = 7.6 (the dashed vertical line denotes to)- Note also that 
the electron's speed decreases slightly in accordance with the previous discussion regarding the 
kinetic energy (see Fig. 2(b)). It should also be noted that this electron deflection time was 
found not at equilibrium conditions but at / = 0.0. If we repeat the measurements (for the same 
electron) during the post-relaxation run (t > 330.0) we obtain t D = 1.0 and /. = 6.8. so the values 
are quite similar to what we had at / = 0.0. In fact, we've found that in general there's no 
appreciable difference between deflection times found at / =0.0 and at / > 330.0. The 
deflection time and the Mean Free Path of a proton were obtained in this case as well; these 
results are listed in Table VII. 

It is interesting to compare these results w ith those given by Spitzer. If we use the 
formula for the deflection time assuming an ideal gas. we get an answer that is too large by 
about an order of magnitude. For a fully ionized gas we must use the equation (Spitzer [3]. p. 
132) 

t D = [8^mt pi (<D(x) - G(x))ln a] ‘ ( 1 8) 

which gives the deflection time for a i ‘test , ‘ particle. Here n is the “field" particles' number 
density (the field particles are all particles other than the ‘Test’' particle), w is the relative speed 
between the test particle and the field particles, p o = (//w’)"'(in this case: ft is the electron- 

proton reduced mass), <P(x) is the error function and G(x) is defined in Spitzer as 

G(x) = (<D(x) - x<&'(x))/2x 2 (19) 
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while x is defined to be 


x = wJfi/2k t T 


( 20 ) 


Spitzer lists values of <t>(x)-G(x) as a function ofx. The quantity InA is a way to take into 

account quantum-mechanical effects. Spitzer tabulates values of InA as a function of 
temperature and density. However, for our problem the density is too high < approx. 9x 1 0 - ' 
electrons/cm 3 ) and so this theory breaks down and he does not list a value of InA appropriate to 
our case. In spite of this, a crude extrapolation of his table allows us to estimate InA - I : 
furthermore, our simulation is purely classical and no quantum effects were taken into account. 
We can take ve ~ (<x % t>+<Vp>)/2 % 5.8 (see Table III below) and from this we get x - 0.6. 
Therefore, according to Spitzer, 0(0.6) - G(0.6) s 0.42 and so plugging into Eq. ( 1 8) we 
obtain t D - 1.3. This is a bit larger than what we obtained, but it’s still in good agreement. 

In Figure 5 w e see a histogram of the instantaneous electron and proton speed 
distributions at the end of this first integration ( t = 330.0). The values of the (instantaneous) 
mean speed and mean squared speed are also shown. The mean speed <v> was obtained first 
and the mean squared speed <v*> was obtained from 

( v : ) = ( v )' +a : (21) 


where o' is the variance of the speed distribution. We represent relevant information about 
these distributions (plus the initial speeds at t =0.0) in the follow ing table: 


Table III: Speed characteristics at t = 330.0 plus initial speeds tail values in atomic units) 


Particle 

Min. speed 

v(t- =0.0) 

<v> 

V< v 2 > 

Max. speed j 

Electron 

1.63 

(12.2322) 

11.33 

12.33 

29.28 

Proton 

0.04 

(0.2854) 

0.26 

OO 

H 

o 

0.70 | 

(el.)/(pr.) 

40.75 

(42.85) 

43.58 

44.03 

41.83 | 


15 



In the last row we divide the electron speed value by the proton speed value. In all cases the 
numbers are very close to the ratio of the square root of their masses \ 1 836.2 =42. 8>- 43. In 
the case of the initial speeds this ratio is exact by definition. Note that even the slowest 
electron is more than twice as fast as the fastest proton. 

Finally in Fig. 6 we plot the positions and velocities of all the particles at the end of 
this first integration. All positions have randomized by this time; note however, the apparent 
formation of a "chain” of particles towards the top of the CB. Again, the velocity vectors of 
most protons are too small to see at this scale, but those of electrons are clearly visible. 

(b) Case N p =1000 particles (Post-Relaxation Run): 

Once the system has relaxed, the program PLSMMD dumps the positions and 
velocities shown in Fig. 6 as initial conditions used in the next program PLSMMD2. The 
reason we have written two separate programs is that once the system is relaxed we search for 
close approaches (an activity that is quite CPU intensive), especially between electrons and 
protons. To speed up the simulation we do not check for conservation of momentum (although 
we still check for energy conservation). 

In Fig. 7 w e see four panels showing the four moments of the electron speed 
distribution. <v>, <v*>, <v J > and <v J > as a function of time ( 0 < / < 440.0 ). The straight 
horizontal lines represent the Maxwell-Boltzmann values and are given by 


< v ) = - 


2k, T 
rm 


(22a) 


, .v 3 k B T 

(»•)- 


m 


(v’H 


r 2klT^ } 


y 7m j 


' ' m~ 


(22b) 

(22c) 

(22d) 


where the q' h moment is given by 
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X 

(v v )= JVg(v)£/v (23) 

0 


and g is the MB speed distribution (Eq. ( 1 ) divided by the number density n, i.e.. g = g / n). T 
= 1 .56x1 0 7 K. is the equilibrium temperature and kg is Boltzmann’s constant. It can be seen 
that at / = 0.0 the moments are far from the MB values and far from relaxation. However, as 
time increases, the moments tend to settle to some "quasi-relaxed” values slightly above the 
MB values. This gives us an estimate of the "relaxation time” for the electrons. We estimate a 
relaxation time for electrons /„ - 84 (2x1 O' 15 s). We note that this value is quite a bit higher 
than that found at other N p 's and this is due to the unusual upward "arching" of the moments at 
the beginning of the integration (especially noticeable in the second moment). The vertical 
dashed line at t = 330.0 separates the "relaxation" run from the "post-relaxation" run, while the 
values in the plots give the mean moments found for the range 330.0 < / < 440.0 (past the 
relaxation time). The follow ing table compares the empirically found values with the MB 
values given by Eqs. (22): 


Table IV: Electron speed distribution moments 


Moment 

Value 

MB (Eqs. 22) 

% Diff. ! 

<v> 

. . 

11.35 

11.22 

*1.1% ; 

<v<> 

151.3 

148.2 

+2.1% j 

i <x?> 

i 

2272 

2217 

-2.5% 

<v 4 > 

\ 

37542 

36610 

*2.5% 


Note that the values are higher by 2% on average than what MB statistics predict. The electron 
speed distribution seems to have reached a ‘‘quasi-equilibrium” temperature that is higher than 
the temperature this problem would suggest. This "electron temperature” can be easily 
obtained from Eq. 22(b), this implies a temperature 2.1% higher than 1 .56x1 0 7 K or 
approximately 1.60xl0 7 K. 

In Fig. 8 w'e see four panels showing the moments of the speed distribution of protons. 
Again, as time increases the moments tend to relax to a “quasi-equilibrium” value, although the 
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time taken is longer than for electrons. It's hard to estimate the equilibrium time from the 
second moment but from third and fourth moment one could estimate t„ ~ 280 (7x1 O' 1 5 s) for 
protons. Just as for electrons, we take the largest value to be the best approximation to the 
relaxation time. i.e.. from the first moment one could estimate in - 190. but the third and fourth 
moments have clearly not vet relaxed. Note that in this case all four moments are lower than 
the MB values (at least after "relaxation"). In Table V we compare the empirical values found 
(again, for the range 330.0 < t < 440.0 ) with the MB values obtained from Eqs. (22): recall 
that the mass of the proton is 1 836.2 au: 


Table V: Proton speed distribution moments 


Moment 

! 

Value 

MB (Eqs. 22) 

% DifT. 

<\> 

0.2594 

0.2618 

-0.9% 

i <r> 

0.0792 

0.0807 

-1.9% 

H 

<r’> 

0.0273 

0.0282 

-3.2% 

| </> 

i 0.0104 

0.0109 

-4.6% | 


Here the empirically found moments are 2.7% lower on average than the MB values. The 
proton temperature implied by <v^> in this case is 1.54xl0 7 K. This is interesting because the 
average of the electron temperature and the proton temperature is very close to the Sun's core 
temperature. It seems that the higher the moment, the further away it is from MB: why this 
should be so is a puzzle. 

What conclusion can w e draw from these plots of the speed distribution? The 
distributions w e have found seem consistent with what Lyman Spitzer wrote in his Physics of 
Fully Ionized Gases , so we will quote him directly (p. 136): 

"We are now in a position to discuss what happens in a proton-electron gas, for 
example, when the velocity distribution is originally arbitrary. We assume that the 
mean kinetic energies of electrons and protons are of the same order of magnitude. 
Collisions of electrons w ith protons will deflect electrons and lead to an isotropic 
velocity distribution, but will not change appreciably the distribution of electron kinetic 
energies. Electron-electron collisions will gradually establish a Maxwellian velocity 
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distribution for the electrons, while proton-proton collisions will yield a corresponding 
velocity distribution for the protons, but at a kinetic temperature that may differ from 
the electron temperature. Finally, equipartition between electrons and protons is 
established by electron-proton collisions."’ 


What we have done is to integrate this "Coulomb” gas forward in time long enough for the 
distribution of both electrons and protons to have relaxed to "quasi-MB“ values, but not so 
long that they have had time to come to thermal equilibrium with each other. Perhaps the 
electrons are "hotter” than protons due to their lighter mass. In any case, according to Spitzer, 
if we had integrated long enough past the "equipartition time” t eq . electrons and protons should 
both have reached the same temperature. The equipartition time is approximately 43 times the 
collision time t c for protons and 1 836 times the collision time for electrons, where the collision 
time in seconds is given by (Spitzer [3]. p. 133): 


_ uaJaF 

c ~ nZ A In A 


(24) 


where A is the particle's mass in terms of the proton mass, n is the number density of a given 
species (per cm 3 ). Z is the charge number (Z = I for both protons and electrons) and from our 
previous discussion. InA = 1 . If we plug in numbers we obtain a collision time for electrons t ce 
- 2xl0 15 s and t cp - 8.6xl0 14 s for protons. Note that the collision time for electrons is 
identical to the relaxation time found previously, while for protons the agreement is off by an 
order of magnitude. Furthermore, note that according to Eq. (24). t ep /t„ = 43 . If we estimate 

the equipartition time to be 1 836 times the electron collision time then we obtain t eq ~ 1 836x80 
= 154224 (in au). With a time step size of l.lxlO 4 this would take 1 .37x1 0 9 time steps (!) to 
achieve, a practical impossibility. 

Finally we are interested in the speed distributions of proton-electron close approaches. 
At a given time step, whenever an electron approached within an arbitrary "small” distance R,„ 
of a proton we recorded the positions and velocities for the pair; we called this a close- 
approach event or simply an “event”. The positions and velocities of the pair were recorded as 
long as the proton-electron distance r y remained smaller than R,„. This was done only during 


19 



the “post-relaxation” run; that is. during 330.0 < / < 440.0 . We define two types of events: 
capture and escape. A “captured" electron has. at a given r . a speed such that it cannot 
escape the Coulomb attraction of the proton while an "escaped” electron has enough speed (at 
a given r ) to escape the Coulomb attraction. The boundary between these two cases is 

obtained from 



(25) 


where p is the electron-proton reduced mass. v cr is the "critical speed”, and e is the softening 
radius; this is just the relative kinetic energy plus the potential energy (corresponding to the 
force function of Eq. (14)). Solving for v cr we obtain 


v 


2 


(26) 


This discussion is of course completely analogous to the concept of parabolic orbits in celestial 
mechanics. At a given electron-proton distance r tJ (with r tj < R in ) we compute the pair's 

relative speed v y = |v ( . - vj : if v,. > v, r then this is an " escape event"; otherw ise it is a "bound 

event". We have found a few cases where an electron was bound in orbit around a proton for 
several time-steps. In such a case, the same electron-proton pair will give rise to several 
consecutive bound-events. Such bound pairs can not last indefinitely because eventually other 
particles" electrical forces will disrupt the bond. In any case, we are mainly interested in 
escape events and the relativ e speed distribution of such events. In Figure 9 we have three 
panels with various speed histograms. We used R m = 5.477 x 10"'" , about 2% of the Wigner- 

Seitz radius. In Figure 9(a) we have a histogram of v* for all close proton-electron 
approaches; there were 6008 “events”. Of these, only 2146 were escape ev ents: in Fig. 9(b) 
we show a histogram of v* for the escape events only. The minimum relative speed squared 

was (v J 2 y ) nun = 345.01 which is v* r at R in = 5.477 xl0‘ 3 while the average relative speed 
squared < v 3 > = 567.05 and the maximum relative speed squared was found to be (v 3 
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= 1 354.0. The histogram “looks" like a MB energy distribution with the zero shifted to the 
right by the amount 345.01 . All of these values are higher than the <v*’> value obtained in 
Table IV ( 1 5 1 .3) because of the proximity of the particles. In other words, the histogram 
shown in Fig. 9(b) is proportional to the energy distribution of passing electrons as "seen" 
from a proton's point of view. For a real Coulomb force and not the modified version we've 
used here, all these values would be even higher. Going back to the alternative nuclear 
reaction e( 3 He. 3 H)v e . recall that we need an electron with at least 18.6 keV = 683.5 hartrees to 
make the reaction feasible. Using the MB value of <v*> from Table IV. the average electron 
energy is approximately //<v>/2 - (0.999)(148.2)/2 = 74.0 hartrees. too low by almost an 
order of magnitude. If we now use the average relative energies of close approaches, w e get 
fj< v 3 >/2 - (0.999)(567.05)/2 = 283.2 hartrees and using the maximum speed we get 

rtv;) max /2 - (0.999)(1354.0) 2 = 676.3. very close to the threshold value. It is clear that if we 

do take into account inter-particle forces, the relative energies increase. Furthermore, had we 
used a "He nucleus instead of a proton in the simulation all these relative energies w ould be 
higher by a factor of Z = 2 due to the higher charge of the 3 He nucleus. The specific v alues 

listed here are dependent on the arbitrarily chosen value of R m = 5.477 x 10~ 3 . A more formal 

analysis would use a value of R,„ as the distance in which the electron and the 3 He nucleus 
"fuse” to become a Tritium nucleus (using the theory of nuclear beta decay). 

If we think of the critical speed as the escape speed of electrons, then we can say that 
electrons w ith v y = v r are on "quasi-parabolic" orbits. Electrons with r v < v. r are on "quasi- 
elliptical" bound orbits (pure, two-bodv orbits in this case are in genera! not closed due to the 
form of the force function, i.e.. Eq. (14)). Meanwhile, electrons with v, ; > v. r are on "quasi- 

hy perbolic” orbits. A particle on a hyperbolic-like orbit would, in the limit of r tJ — ► x. have a 
certain “excess speed" above zero (if it were zero we would have a parabolic-like orbit). For 
every event shown in Fig. 9(b) we have subtracted an amount corresponding to a parabolic 
orbit (given by Eq. (26)) and show the results in the last histogram. Fig. 9(c): this is a 
histogram of the “excess speed squared”. The figure is a half-Gaussian with a standard 
deviation of cr= 1 13 and a mean value of < vj> = 1 15.17. It was expected that perhaps the 

values of < v 3 > for the excess speed would be higher than the mean squared speed for an ideal 
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gas (i.e.. < v* > = 148.2), but this was not the case. We checked to see if this deficiency could 
be explained by Debye screening, but screening is a very small effect at this close range 
( r. < R,„ ) and does not account for the smallness of < \’l >. 

An attempt was also made to obtain histograms for electron-electron and proton-proton 
close approaches but these happened too infrequently to obtain useful data. 

(c) Summary of Results: 

In this section we summarize our results using S’ p = 512, 1000, 1728 particles, which 
are presented in the form of tables. Again, the time step size is Ai = l.lxl 0 -4 , and the outer 
"cut-off' radius R c = 0.67. 

In Table VI we show the electron deflection times plus the Mean Free Paths. The 
analytical mean free path is obtained from /. = <v> S(B to =(11 -22X 1.3)= 1 4.6. An asterisk 
means that we are taking the average value obtained starting from t= 0.0 and from t = 330.0. 
since, as we have seen for the N p = 1 000 case, the values do not vary much. The values are 
somewhat consistent w ith each other. It should also be noted that in all cases the Mean Free 

Paths are larger than the sides of the largest CB's side (Z,^ = r/(l728) ’ s 4.0). 


Table VI: Electron Deflection Times & Mean Free Paths 


N p 

Deflection Time t D 

Mean Free Path k 

512* 

1.9 

26.6 

1000* 

0.9 

7.2 

1728* 

1.1 

9.6 

Spitzer (Eq. 18) 

1.3 

14.6 


In Table VII we show the deflection times plus the Mean Free Paths for protons. These 
values were measured only at / = 0.0. Unfortunately we could not obtain values for N p = 1 728. 
However, it is apparent that the measured deflection times are higher than Spitzer s 




predictions. Theoretically we expect the Mean Free Paths of electrons and protons to be the 


same. 


Table VII: Proton Deflection Times & Mean Free Paths 


.% 

Deflection Time t D 

Mean Free Path A 

512 

103.0 

i6.i ; 

1000 

206.0 

52.1 | 

1728 

na 

na 

Spitzer (Eq. 18) 

55.7 

14.6 


In Table VIII we show the values of the estimated relaxation times t n for protons and 
electrons. These relaxation times are a measure of how long it takes the speed distribution to 
relax from FCC lattice ty pe initial positions and delta functions as the initial speed 
distributions. The means are also computed: as expected, the relaxation time for electrons is 
much shorter than that for protons due to the latter's higher mass and lower mean speeds. 


Table VIII: Relaxation Times 


% 

Electrons 

Protons 

512 

21 

150 

1000 

84 

280 

1728 

14 

300 

Mean 

39.7 

243.3 


In Table IX we present the observed values of the first four moments of the electron 
speed distribution. Recall that these values were obtained only for the range 
330.0 < / < 440.0 , i.e., the “post-relaxation'’ run (330 < t < 424.28 for the case \ p = 1 728). 
Mean values were obtained and percent differences between these mean values and the MB 
values are also shown. It seems that the higher the number of particles, the higher the value of 
the moment. However, if we make the assumption that the higher the number of particles the 
more “accurate” the simulation is, the values seem to be converging “up" to some asy mptote as 
the number of particles get higher. This is shown in Figure 10 for the second moment: the 




other moments show this type of behav ior as well. It w ould be hard to justify a fit to find this 
convergent value with only three points, however 


Table IX: Electron speed distribution moments (330.0 <t < 440.0 ) 


% 

<v> 

<v r ’> j <v*> 

<v 4 > 

512 

11.31 

149.70 

2231.00 

36640.00 

1000 j 11.35 

151.30 

2272.00 

37540.00 

1728* 

11.38 

151.89 

2283.40 

37724.00 

Mean 

1 1 .35 

150.96 

2262.13 

37301.33 

MB 

11.22 

148.20 

2217.00 

36610.00 

%Difference 

1.1% 

1.9% 

2.0% 

i 

1.9% 


An asterisk means that the time range is 330 < t < 424.28. 


In Table X w e show the values of the first four moments of the proton speed 
distribution for 330.0 < t < 440.0 (330 < t < 424.28 for the case \ P = 1 728). Again we take 
the average and compare it to the MB value. Here, the average values are lower, in contrast to 
the electron case. Interestingly, the higher the moment the further away the values are from 
MB. In Figure II we show the values of the second moment as a function of the number of 
particles. Again, the empirical values seem to be converging "down" to some asymptotic 
\alue; the other moments seem to converge also. 

Finally in Table XI we show the "mean relative speeds squared", raw (that is. just what 
was directly measured <v' > (raw) ) and the "mean relative excess speed squared" <v* > (eycess) . 
The relationship between these two is 

^ V y >( ivcess ) — <V V >( raw) — (27) 

where v cr depends on the relative electron-proton distance r tJ and is given by Eq. (26). Note 

that the mean squared excess speed raw (1 13.66). is smaller than the mean squared speed 
assuming a MB distribution (148.20). 
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Table X: Proton speed distribution moments (330.0 < t < 440.0 ) 


V 

<v> 

, 

<r> 

</> 

<v 4 > 

512 

0.2600 

0.0799 

0.0277 

0.0105 

1000 

0.2594 

0.0792 

0.0273 

0.0104 

1728* 

0.2601 

0.0791 

0.027 

0.01 

Mean 

0.2598 

0.0794 

0.0273 

0.0103 

MB 

0.2618 

0.0807 

0.0282 

0.0109 

%DifT 

-0.8% 

-1.6% 

-3.1% 

-5.5% 


An asterisk means that the time range is 330 < t < 424.28. 


Table XI: Mean relative speeds squared, raw and excess 


( R m = 5.477 x 1 0'\ 330.0 < t < 440.0 ) 


N, 

< x \j > ( raw) 

< v .~ > (excas) 

512 

570.50 

113.50 

1000 

567.05 

115.17 

1728* 

566.44 

112.31 

Mean 

568.0 

113.66 


An asterisk means that the time range is 330 < t < 424.28. 


V. CONCLUSIONS/SUMMARY 

We have integrated a simplified model of a plasma forward in time using a Molecular 
Dynamics approach. The total number of particles was varied and results with different \ p 
w ere compared. This plasma is a fully ionized gas consisting of protons and free electrons and 
has the same temperature and density as is found in the core of the Sun. We wished to 
investigate whether inter-particle forces (the Coulomb force) in this multi-body problem might 
be a source of deviations from MB statistics predictions, which are strictly true for an ideal gas. 
Specifically we wished to see if the reaction given by Eq. (4), e'( 3 He, 3 H)v e , [Q = -1 8.6 keV] 
might be made more likely than standard solar theory predicts due to multi-body effects. 

The plasma was integrated from non-equilibrium initial conditions (FCC lattice for positions 
and a delta function for the speed distribution) until quasi-equilibrium was reached. During 
this time we checked the deflection times, the mean free paths and the relaxation times (from 
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the first four moments) for electrons and protons. The results were somewhat consistent w ith 
analytical results as provided by Spitzer [3], especially for electrons. 

After the system was completely "quasi-relaxed", we again propagated the system 
forward in time and collected statistics on the first four moments of the proton and electron 
speed distributions as well as statistics on close inter-particle approaches. The first four 
moments of the speed distribution were recorded for the duration of the post-relaxation run and 
the average results were compared to the MB values. A curious phenomenon is that, although 
we started protons and electrons with the same kinetic energies, electrons quickly gained 
kinetic energy at the expense of protons. Therefore, electrons were slightly "hotter” (+ 1 .9%) 
and protons were slightly cooler (-1.6%) than the system temperature of 1.56xl0 7 K (see 
Tables IX and X). However, this is not at all inconsistent with the behavior of a fully ionized 
gas. as shown by Spitzer [3], Note that the average of the electron ( 1 .59x I O' K) and proton 
temperatures (1.53x!0 7 K) is very close to the system temperature. To have the plasma in true 
equilibrium we would have to have integrated the system past the "equipartition time”, where 
protons and electrons would have the same temperature. As we have seen, this would be out of 
the question w ith the present approach. This is w hy we prefer to call the plasma "quasi- 
relaxed” and not relaxed at the end of the integration. 

A way we used to measure deviation from MB was to look at the first four moments 
during "quasi-equilibrium (t > 330.0) and compare the values with MB. The electron moments 
were higher (+1.7% on average) than MB and the proton moments were lower (-2.8% on 
average) than MB. However, it is likely that the observed deviations from MB are consistent 
with the plasma not being in complete thermal equilibrium. 

What about close approaches between electrons and protons? We do have to use the 
total relative energy and not just the kinetic energies. For example in the N p = 1000 particles 
case, we had that the mean relative kinetic energy was <k> = 283.2 hartrees, but the mean 

relative potential energy in this case was <u> = - Jr* + ~e 2 J” = -226.4 hartrees for a total 

mean relative energy of <e> = 56.8 hatrees, lower than even the MB relative energy of 74.0 
hartrees and quite low compared to the threshold energy- for the reaction e( 3 He, 3 H)v e , [Q = - 
683.5 hartrees]. Note that, relative to the MB case, the inter-particle forces seem to have 


26 



actually decreased the mean relative energies, possibly making the rate of the reaction 
e'(’He.'H)v e less likely than in the pure MB case. As previously stated, the mean squared 
excess speed raw (1 13.66). is smaller than the MB mean squared speed (148.20). 

We conclude that the observed deviations from a MB distribution are probably not 
significant enough to make the alternative reaction e ( ? He. 3 H)v e more likely than in a pure 
MB case. Indeed, the observed deviations are in accordance with the expected behavior of a 
plasma, as described by Spitzer. Deviations from MB due to electron degeneracy (protons 
are practically non-degenerate) were also explored analytically as an aside project. At the 
temperature and pressure found at the core of the Sun. the mean electron energy is about 3% 
higher in the degenerate case than in the MB case. Again, this is probably too small an effect 
to make a difference. It seems that, in light of recent inv estigations regarding neutrino 
"(lavor changes’", a possibly better explanation for the deficiency of neutrinos from the Sun 
may be neutrino oscillations, or the MSW effect [ 1 7, 1 8]. Very briefly, neutrinos exist in 
three types: electron-neutrino ( v,). which is the only type that is produced in the Sun, the 
muon-neutrino ( v M ) and the tau-neutrino ( vy). Recent experiments at the Super-Kamiokande 
neutrino detector in Japan have suggested that oscillations indeed take place in the case of the 
tau and muon neutrinos and possibly in the case of electron-neutrinos [19]. This implies that 
neutrinos have mass, contrary to standard particle physics, which assumes zero mass for 
them. If oscillations do take place, as (electron) neutrinos from the Sun travel to Earth, they 
could sw itch to either of the other two types (tau and muon). With the exception of the 
Sudbury Neutrino Observatory in Canada, which can detect all three flavors and w ill test the 
MSW effect, present-day “neutrino telescopes” can only detect one or two types of neutrinos 
at a time. If neutrinos oscillate, the deficiency of solar neutrinos may be only an apparent 
one, as the detectors only detect neutrinos of one flavor (electron) or the other, thus possibly 
explaining the solar neutrino problem. Therefore, it may be that the answer to the solar 
neutrino problem may lie in particle physics, not statistical mechanics. 
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VI. FIGURES AND FIGURE CAPTIONS 


FIGURE 1: Initial conditions (positions and velocities) for the particular case N p = 1000 
particles. Initial positions correspond to an FCC lattice. Only the velocity vectors of the 
electrons can be seen. 

FIGURE 2: Behav ior of the system energy for the case S p = 1000. (a) System energy 
(kinetic, total and potential), (b) System energy but with the potential energy increased by 
72000 hartrees. (c) System energy separated into proton and electron components, (d) 
Percent error in the total energy. 

FIGURE 3 : Behavior of the system momentum for the case S p = 1 000. (a) Total and x, y, r 
components of the system momentum, (b) Electron x. y. z momentum components. (c) 
Proton x, y, r momentum components. 

FIGURE 4: Behavior of a given particle (in this case particle number 18, an electron) as a 
function of time for the case N p = 1 000. (a) Position as a function of time (shown in the form 
of components), (b) Velocity as a function of time (shown in the form of components). (c) 
Cosine of the angle betweeen the initial velocity' vector v 0 and the velocity vector as a 
function of time v(t). Circle represents the deflection time (here 0.80). (d) Speed of the 
particle as a function of time. v(t) vs. t. The dashed vertical line corresponds to the deflection 
time and the integral of vfi) from zero to the deflection time gives the mean free path (here 
equal to 7.6 time units). 

FIGURE 5: Speed histograms for electrons (panel a) and protons (panel b) for the case N p = 
1000 at the end of the “relaxation run" (/ = 330.0). 

FIGURE 6: Final conditions (positions and velocities) for the particular case N p = 1000 
particles at the end of the “relaxation run” (t = 330.0). Only the velocity vectors of the 
electrons can be seen. 

FIGURE 7: First four moments of the electron distribution as a function of time for the case 
N p = 1000. The horizontal lines represent the Maxwell-Boltzmann values. 
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FIGURE 8: First four moments of the proton distribution as a function of time for the case 
X p = 1000. The horizontal lines represent the Maxwell-Boltzmann values. 

FIGURE 9: Distributions of squared speed for close electron-proton approaches X p = 1000 
and R,„ = 0.005477. These “events'* cover the period of the “post-relaxation" run 330.0 < t < 
440.0. (a) Histogram of the squared speed for all events, (b) Histogram of the squared 
speed for escape events only, (c) Histogram of the squared speed minus the critical speed. 

In other words, histogram of the excess speed squared. 

FIGURE 10: Empirical values of the mean squared speed <v'> for electrons as a function of 
the number of particles (cases X p = 512. 1 000. 1 728). It can be seen that the mean value is 
higher than the Maxwell-Boltzmann value. 

FIGURE 1 1 : Empirical values of the mean squared speed <v’> for protons as a function of 
the number of particles (cases X p = 512. 1 000. 1 728). It can be seen that the mean value is 
lower than the Maxwell-Boltzmann value. 
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APPENDIX: 

Inputs to and programs PLSMMD.FOR & PLSMMD2.FOR 

The first thing to decide before running these two programs is the number of 
panicles. Recall that the number of particles is restricted by the condition N p = 8 i 3 where / = 
1 .2,3 .4.5.6.... etc. or A r p = 8, 64, 216, 512, 1000, 1728. In this “experiment” we only 
considered 5 12 < N p < 1 728. Once the number of particles has been decided upon, the 
FORTRAN include file PLASMA.INC is edited w ith a text editor and the number of 
particles is entered, where it is defined a PARAMETER: 

T ....IRE TOTAL NUMBER OF PARTICLES 

- NP=64, 216, 512, 1000, 1728, 8*I~3 (I=l,2,3...) 

INTEGER NP 
PARAMETER (NP = 512) 

The program PLSMMD.FOR is then compiled. Now we have the executable, which needs 
an input file stating the number of time-steps, the time-step size. etc. The contents of this 
input file LN.DAT are as follows: 


0.26 

The Wigner-Seitz radius, in atomic units (au) 

1.56e7 

The temperature, in Kelvin 

1000000 

The number of time steps 

1.1 e-4 

The time-step size (au) 

0.67 

"Outer*' cut-off radius (au) 

18 

Particle number to follow (output in file TXYZ.DAT) 

100 

How often to WTite energies file ENERGY.DAT (# of time-steps) 



15 

0.5 


How often to update the neighbor lists (# of time-steps) 
Max. % error allowed in energy 


Once the program PLSMMD.FOR is finished, it outputs the file IN2.DAT. which is the 
input file for the program PLSMMD2.FOR. The contents of the file 1N2.DAT are as 
follows: 


512 

330.0 
0.26 

15600000.000000 
1000000 
1.10D-04 

0.67 

18 

100 

15 

1 

0.00547722557052 

36959.274946828 

-290.94705096300 

0.50000000000000 


Total # of particles (used as a check) 

Time at the end of the first run [computed by PLSMMD.FOR] 
Wigner-Seitz radius, in atomic units (au) 

Temperature, in Kelvin 
Total # of time-steps 
Time-step size (au) 

"Outer" cut-ofT radius (au) 

Particle number to follow (output in file TXYZ2.DAT) 

How often to write energies file ENERGY2.DAT (# of time-steps) 
How often to update the neighbor lists (# of time-steps) 

How often to check close approaches (# of time-steps) 

How close panicles can approach to record relative pos., vel. 

Total energy (hanree) [computed by PLSMMD.FOR] 

Corr. to the potential energy (hartree) [computed by PLSMMD.FOR] 
Max. % error allowed in energy' 


The first number is used as an internal check to make sure the two programs are compiled to 
the same number of panicles. Note that some of the numbers are computed by the program 
PLSMMD.FOR and the user is not allowed to change them (for example, the time at the end 
of the PLSMMD.FOR run), lest there be inconsistencies in the two programs. 



